function [suppr,freq] = twoshiftpd(N,D1,D2,M,wr,sc)
% twoshiftpd(N,D1,D2,M,wr,sc)
% 	Parameters:
% 		N:	number of data points to include in time series
% 		D1:	first number of seconds to displace time series over
% 		D2:	second number of seconds to displace time series over
% 		M:	kernel size
% 		wr:	type of noise to generate ('w' = white noise, 'r' = 1/f noise, 'b' = f noise)
% 		sc:	coefficient to scale noise
% 	Outputs:
% 		suppr	y-values for a suppression plot
% 		freq 	x-values for a suppression plot	
% 	Purpose of this code:  This code will produce a time series of noise in the time domain and will pass it through a lowpass filter then will shift 
% 	it once D1 seconds then again for D2 seconds by interpolation.  Then this code will take the same time series and shift it in all one shift again by 
% 	interpolation.  The difference between the two timeseries will be found.  The power spectrum of the result will be divided by the power spectrum
% 	of the original timeseries and a square root taken to arrive at the y-values for a suppression plot.  This plot is provided by this code.


%hold off;
%% generate sequence of white noise to be sampled
%pool = whitenoise(PS.*PRes+1);
%% sample above sequence

%ts(1) = pool(1);
%for j = 2:(PS.*rate)
%    ts(j) = pool((j.*PRes+1)./rate);
%end
%length(ts)
ts = whitenoise(N,wr,sc);
ts = lowpass(ts);
[psorig freq] = mypsd(ts);

tsonce = FDtest(ts,(D1+D2),M);
psa = mypsd(tsonce);
%% shift timeseries by S1 and then by S2
tsd1 = FDtest(ts,D1,M);
tsd2 = FDtest(tsd1,D2,M);

ps = mypsd(tsd2);

%% First Plot
%loglog(1:length(psorig),psorig,'b');
%hold on;
%err = abs(ps - psa)./abs(psa);
%loglog(1:length(ps),err,'g');
err2 = abs(tsd2-tsonce);
errps = mypsd(err2);
%figure
%loglog(freq,psorig,freq,errps);
%figure
%plot(1:N,ts,1:N,err2);
%plot(1:N,tsd2,1:N,tsonce,1:N,err2);

suppr = sqrt(psorig./errps); % taking sqrt to go from power to amplitude
%figure
%loglog(freq,suppr);

xnn = abs(ps - psa);
supp = psa./xnn;
%loglog(1:length(ps),supp,'r');
